Group 1

rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Group == 1])
#rnaseq.subset <- subset.PDX(rnaseq.subset, rnaseq.subset$meta[,Sample.name!="HCI011-T1-P01-18.9413R1-T1"])
treatment <- rnaseq.subset$meta[,factor(Treatment.arm,levels=c("vehicle","GDC0032"))]
d <- model.matrix(~ 0 + treatment)
d

  treatmentvehicle treatmentGDC0032
1                1                0
2                1                0
3                1                0
4                1                0
5                0                1
6                0                1
7                0                1
8                0                1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment
[1] "contr.treatment"

cont.matrix <- matrix(c(-1,1), dim(d)[2], 1)
colnames(cont.matrix) <- c("GDC.minus.control")
rownames(cont.matrix) <- colnames(d)
cont.matrix

                 GDC.minus.control
treatmentvehicle                -1
treatmentGDC0032                 1

#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,cont.matrix = cont.matrix)

#Top DE genes table
DE <- topTable.annotated.DE(DE.data$efit,p.val = 1)

'select()' returned 1:many mapping between keys and columns

if(nrow(DE[adj.P.Val<0.05])!=0){
  DE <- DE[adj.P.Val<0.05]
  table(DE[,direction])
  DE[direction=="up.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
  DE[direction=="down.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
#Pheatmap w/ batch correction
  i <- match(DE$ENSEMBL,rownames(DE.data$v$E))
  y <- removeBatchEffect(DE.data$v$E[i,],design = model.matrix(~drug), batch= pdx)
  labels.col <- rep("",nrow(rnaseq$meta))
  labels.row <- DE$SYMBOL
  ann.col <- as.data.frame(rnaseq$meta[,.(pdx,drug)])
  rownames(ann.col) <- rnaseq$meta[,Sample.name]
  pheatmap(y,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
} else {
  head(DE,10)
}

            ENSEMBL      logFC  AveExpr         t      P.Value adj.P.Val         B direction ENTREZID    SYMBOL                                              GENENAME
 1: ENSG00000204388  1.0760352 9.253658  4.832736 0.0006147231 0.3272351 -1.519333    up.reg     3304    HSPA1B         heat shock protein family A (Hsp70) member 1B
 2: ENSG00000236675 -0.8133502 5.721730 -4.733665 0.0007165805 0.3272351 -1.574151  down.reg     <NA>      <NA>                                                  <NA>
 3: ENSG00000183011 -0.8991874 6.269441 -4.568031 0.0009288204 0.3272351 -1.620282  down.reg    84316     NAA38 N(alpha)-acetyltransferase 38, NatC auxiliary subunit
 4: ENSG00000145494 -0.7980749 6.345503 -4.512080 0.0010147699 0.3272351 -1.652281  down.reg     4726    NDUFS6             NADH:ubiquinone oxidoreductase subunit S6
 5: ENSG00000237991 -1.0496058 6.159797 -4.452303 0.0011159358 0.3272351 -1.719579  down.reg     <NA>      <NA>                                                  <NA>
 6: ENSG00000151208  1.0660589 5.918073  4.490701 0.0010497920 0.3272351 -1.728342    up.reg     9231      DLG5                  discs large MAGUK scaffold protein 5
 7: ENSG00000077809  0.8711585 6.583563  4.340502 0.0013347542 0.3272351 -1.771400    up.reg     <NA>      <NA>                                                  <NA>
 8: ENSG00000105669 -0.7967949 7.069990 -4.264799 0.0015082447 0.3272351 -1.807395  down.reg    11316      COPE              coatomer protein complex subunit epsilon
 9: ENSG00000125356 -0.6583621 6.864270 -4.197563 0.0016822276 0.3272351 -1.863827  down.reg     4694    NDUFA1             NADH:ubiquinone oxidoreductase subunit A1
10: ENSG00000180389 -0.7912983 6.400057 -4.216161 0.0016320902 0.3272351 -1.873863  down.reg   432369 ATP5F1EP2          ATP synthase F1 subunit epsilon pseudogene 2

#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Treatment.arm) + geom_boxplot() + stat_compare_means(aes(label=paste0("p.adj = ",..p.adj..)),method="t.test") + theme_tufte() + theme(legend.position = "top")

## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Treatment.arm)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")

#GSEA
## C6
gsea <- GSEA.run(DE.data$v,d,cont.matrix[,1],Hs.c6,0.05)

'select()' returned 1:many mapping between keys and columns

gsea[Direction=="Up",.(rn,NGenes,Direction,FDR)]

                                  rn NGenes Direction          FDR
 1:                    TGFB_UP.V1_UP    133        Up 5.785652e-05
 2:                    STK33_NOMO_UP    210        Up 2.290397e-04
 3:                         STK33_UP    201        Up 3.175661e-04
 4:                    LTE2_UP.V1_DN    174        Up 3.972324e-04
 5:                     STK33_SKM_UP    192        Up 9.710797e-04
 6:               BCAT_BILD_ET_AL_DN     38        Up 9.766632e-04
 7:                         EIF4E_DN     90        Up 2.000560e-03
 8:                    BMI1_DN.V1_UP    101        Up 2.347305e-03
 9:           BMI1_DN_MEL18_DN.V1_UP     95        Up 2.444336e-03
10:                       TBK1.DF_DN    266        Up 5.350087e-03
11:                    CAMP_UP.V1_DN    155        Up 5.671311e-03
12:               CSR_EARLY_UP.V1_DN    121        Up 6.646370e-03
13:                    PIGF_UP.V1_UP    152        Up 8.964812e-03
14:                    EGFR_UP.V1_UP    168        Up 8.964812e-03
15:                     MEK_UP.V1_DN    172        Up 1.038029e-02
16:                     RAF_UP.V1_DN    177        Up 1.038029e-02
17:                   MEL18_DN.V1_UP     96        Up 1.299166e-02
18:                   RPS14_DN.V1_UP     97        Up 1.354666e-02
19:               PRC2_EZH2_UP.V1_DN    114        Up 1.940560e-02
20:                    PDGF_UP.V1_DN     70        Up 1.940560e-02
21:                PDGF_ERK_DN.V1_DN    130        Up 2.272031e-02
22:                     RAF_UP.V1_UP    167        Up 2.272031e-02
23:                 RB_P130_DN.V1_DN    121        Up 2.272031e-02
24:                    KRAS.DF.V1_DN    129        Up 2.272031e-02
25:             HINATA_NFKB_IMMU_INF      7        Up 2.297676e-02
26:                 MTOR_UP.N4.V1_DN    146        Up 2.652461e-02
27:                    PDGF_UP.V1_UP    127        Up 3.046860e-02
28:                    E2F1_UP.V1_DN    148        Up 3.196333e-02
29:                   HOXA9_DN.V1_UP    148        Up 3.495692e-02
30:                PRC1_BMI_UP.V1_DN     86        Up 3.584935e-02
31:                PRC2_EED_UP.V1_UP    108        Up 3.584935e-02
32:                     ALK_DN.V1_DN     64        Up 4.019820e-02
33:                    RELA_DN.V1_UP     92        Up 4.672090e-02
34:                     MEK_UP.V1_UP    170        Up 4.672090e-02
35:                    IL21_UP.V1_DN     74        Up 4.672090e-02
36:                    ERB2_UP.V1_DN    185        Up 4.672090e-02
37:                     MYC_UP.V1_DN    116        Up 4.672090e-02
38: SINGH_KRAS_DEPENDENCY_SIGNATURE_     18        Up 4.672090e-02
                                  rn NGenes Direction          FDR

gsea[Direction=="Down",.(rn,NGenes,Direction,FDR)]

                  rn NGenes Direction          FDR
1:          EIF4E_UP     77      Down 0.0002290397
2: CSR_LATE_UP.V1_UP    142      Down 0.0002463503
3:  RB_P107_DN.V1_UP    116      Down 0.0002496484
4:  RB_P130_DN.V1_UP    100      Down 0.0194056047
5:      MYC_UP.V1_UP    142      Down 0.0227203071
6:    RPS14_DN.V1_DN    155      Down 0.0236125612

Group 2

rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Group == 2])
treatment <- rnaseq.subset$meta[,factor(Treatment.arm,levels=c("vehicle","GDC0032"))]
d <- model.matrix(~ 0 + treatment)
d

  treatmentvehicle treatmentGDC0032
1                1                0
2                1                0
3                1                0
4                0                1
5                0                1
6                0                1
7                0                1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment
[1] "contr.treatment"

cont.matrix <- matrix(c(-1,1), dim(d)[2], 1)
colnames(cont.matrix) <- c("GDC.minus.control")
rownames(cont.matrix) <- colnames(d)
cont.matrix

                 GDC.minus.control
treatmentvehicle                -1
treatmentGDC0032                 1

#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,3,cont.matrix = cont.matrix)

#Top DE genes table
DE <- topTable.annotated.DE(DE.data$efit,p.val = 1)

'select()' returned 1:many mapping between keys and columns

if(nrow(DE[adj.P.Val<0.05])!=0){
  DE <- DE[adj.P.Val<0.05]
  table(DE[,direction])
  DE[direction=="up.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
  DE[direction=="down.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
#Pheatmap w/ batch correction
  i <- match(DE$ENSEMBL,rownames(DE.data$v$E))
  y <- removeBatchEffect(DE.data$v$E[i,],design = model.matrix(~drug), batch= pdx)
  labels.col <- rep("",nrow(rnaseq$meta))
  labels.row <- DE$SYMBOL
  ann.col <- as.data.frame(rnaseq$meta[,.(pdx,drug)])
  rownames(ann.col) <- rnaseq$meta[,Sample.name]
  pheatmap(y,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
} else {
  print("No DE genes found")
  head(DE,10)
}

[1] "No DE genes found"
            ENSEMBL      logFC  AveExpr         t      P.Value adj.P.Val         B direction ENTREZID  SYMBOL                                                                        GENENAME
 1: ENSG00000174021  1.0508825 6.847082  6.713034 7.674007e-05 0.6224003 -3.671849    up.reg     2787    GNG5                                                       G protein subunit gamma 5
 2: ENSG00000164128 -1.1130122 5.649225 -5.955737 1.919213e-04 0.9268516 -3.790498  down.reg     4886   NPY1R                                                      neuropeptide Y receptor Y1
 3: ENSG00000204569  0.8046126 6.989129  4.701910 1.037803e-03 0.9268516 -3.835360    up.reg     5514 PPP1R10                                     protein phosphatase 1 regulatory subunit 10
 4: ENSG00000164442 -0.9972152 5.227045 -5.581024 3.105660e-04 0.9268516 -3.852576  down.reg    10370  CITED2 Cbp/p300 interacting transactivator with Glu/Asp rich carboxy-terminal domain 2
 5: ENSG00000186918 -0.6304149 5.900330 -4.502097 1.386601e-03 0.9268516 -3.907247  down.reg    55893  ZNF395                                                         zinc finger protein 395
 6: ENSG00000104267 -0.8459846 6.814805 -4.081274 2.600769e-03 0.9268516 -3.919892  down.reg      760     CA2                                                            carbonic anhydrase 2
 7: ENSG00000144802 -1.1270393 5.645451 -4.275787 1.938616e-03 0.9268516 -3.952431  down.reg    64332  NFKBIZ                                                             NFKB inhibitor zeta
 8: ENSG00000074800 -0.4960966 8.821536 -3.652362 5.062147e-03 0.9268516 -3.974104  down.reg     2023    ENO1                                                                       enolase 1
 9: ENSG00000135127  0.8043419 6.155308  3.951610 3.172693e-03 0.9268516 -3.976310    up.reg    92558  BICDL1                                                BICD family like cargo adaptor 1
10: ENSG00000104321 -1.8460421 7.281383 -3.667565 4.942057e-03 0.9268516 -3.983110  down.reg     8989   TRPA1                transient receptor potential cation channel subfamily A member 1

#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Treatment.arm) + geom_boxplot() + stat_compare_means(aes(label=paste0("p.adj = ",..p.adj..)),method="t.test") + theme_tufte() + theme(legend.position = "top")

## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Treatment.arm)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")

#GSEA - c6
gsea <- GSEA.run(DE.data$v,d,cont.matrix[,1],Hs.c6,0.05)

'select()' returned 1:many mapping between keys and columns

gsea[Direction=="Up",.(rn,NGenes,Direction,FDR)]

                                  rn NGenes Direction          FDR
 1:          ESC_V6.5_UP_EARLY.V1_DN    116        Up 8.734097e-06
 2:                    BMI1_DN.V1_UP    101        Up 1.146764e-03
 3:           BMI1_DN_MEL18_DN.V1_UP     94        Up 1.299194e-03
 4:                    TGFB_UP.V1_UP    133        Up 3.186118e-03
 5:                KRAS.600_UP.V1_DN    134        Up 3.186118e-03
 6:                    PTEN_DN.V2_UP     78        Up 4.432734e-03
 7:           KRAS.AMP.LUNG_UP.V1_DN     74        Up 5.805192e-03
 8:                KRAS.300_UP.V1_DN     72        Up 8.586000e-03
 9:               KRAS.LUNG_UP.V1_DN     68        Up 1.156810e-02
10:                   MEL18_DN.V1_UP     95        Up 1.211232e-02
11:                    SNF5_DN.V1_DN     92        Up 3.277517e-02
12:             KRAS.KIDNEY_UP.V1_DN     60        Up 3.277517e-02
13:                 CAHOY_ASTROGLIAL     59        Up 3.277517e-02
14:                    RELA_DN.V1_DN     58        Up 3.277517e-02
15:                     MEK_UP.V1_UP    171        Up 3.277517e-02
16:               CYCLIN_D1_UP.V1_UP    127        Up 3.874861e-02
17:           CAHOY_OLIGODENDROCUTIC     69        Up 4.393997e-02
18:                     P53_DN.V1_UP    147        Up 4.393997e-02
19: SINGH_KRAS_DEPENDENCY_SIGNATURE_     18        Up 4.393997e-02
20:              RAPA_EARLY_UP.V1_UP    115        Up 4.393997e-02
21:                     MYC_UP.V1_DN    115        Up 4.589087e-02
22:                  VEGF_A_UP.V1_UP    129        Up 4.589087e-02
23:                     AKT_UP.V1_DN    132        Up 4.589087e-02
                                  rn NGenes Direction          FDR

gsea[Direction=="Down",.(rn,NGenes,Direction,FDR)]

                rn NGenes Direction        FDR
1:  HOXA9_DN.V1_DN    160      Down 0.03277517
2:      TBK1.DF_DN    268      Down 0.03831569
3: VEGF_A_UP.V1_DN    167      Down 0.04466353

Group 3

rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Group == 3])
treatment <- rnaseq.subset$meta[,factor(Treatment.arm,levels=c("vehicle","GDC0032"))]
d <- model.matrix(~ 0 + treatment)
d

   treatmentvehicle treatmentGDC0032
1                 1                0
2                 1                0
3                 1                0
4                 1                0
5                 1                0
6                 0                1
7                 0                1
8                 0                1
9                 0                1
10                0                1
11                0                1
12                0                1
13                0                1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment
[1] "contr.treatment"

cont.matrix <- matrix(c(-1,1), dim(d)[2], 1)
colnames(cont.matrix) <- c("GDC.minus.control")
rownames(cont.matrix) <- colnames(d)
cont.matrix

                 GDC.minus.control
treatmentvehicle                -1
treatmentGDC0032                 1

#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,5,cont.matrix = cont.matrix)

#Top DE genes table
DE <- topTable.annotated.DE(DE.data$efit,p.val = 1)

'select()' returned 1:many mapping between keys and columns

if(nrow(DE[adj.P.Val<0.05])!=0){
  DE <- DE[adj.P.Val<0.05]
  table(DE[,direction])
  DE[direction=="up.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
  DE[direction=="down.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
#Pheatmap w/ batch correction
  i <- match(DE$ENSEMBL,rownames(DE.data$v$E))
  y <- removeBatchEffect(DE.data$v$E[i,],design = model.matrix(~drug), batch= pdx)
  labels.col <- rep("",nrow(rnaseq$meta))
  labels.row <- DE$SYMBOL
  ann.col <- as.data.frame(rnaseq$meta[,.(pdx,drug)])
  rownames(ann.col) <- rnaseq$meta[,Sample.name]
  pheatmap(y,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
} else {
  print("No DE genes found")
  head(DE,10)
}

[1] "No DE genes found"
            ENSEMBL      logFC  AveExpr         t      P.Value adj.P.Val           B direction ENTREZID SYMBOL                                                   GENENAME
 1: ENSG00000163995  0.6837759 3.410971  5.496805 3.477881e-05 0.5074877  0.17215334    up.reg    84448 ABLIM2                  actin binding LIM protein family member 2
 2: ENSG00000111300  0.4156489 5.496378  4.531102 2.726366e-04 0.5074877  0.10500080    up.reg    80018  NAA25      N(alpha)-acetyltransferase 25, NatB auxiliary subunit
 3: ENSG00000102753  0.3786802 6.106180  4.424698 3.438274e-04 0.5074877  0.03975911    up.reg     3839  KPNA3                                karyopherin subunit alpha 3
 4: ENSG00000253729  0.4086519 6.749566  4.330204 4.227513e-04 0.5074877 -0.05735691    up.reg     5591  PRKDC       protein kinase, DNA-activated, catalytic polypeptide
 5: ENSG00000204389 -0.5307251 7.717875 -4.208828 5.516668e-04 0.5074877 -0.22188442  down.reg     3303 HSPA1A              heat shock protein family A (Hsp70) member 1A
 6: ENSG00000139055 -0.4244053 5.203128 -4.308408 4.434221e-04 0.5074877 -0.26761439  down.reg   121506  ERP27                           endoplasmic reticulum protein 27
 7: ENSG00000204922 -0.4767766 4.225305 -4.500622 2.913477e-04 0.5074877 -0.27855530  down.reg   790955  UQCC3 ubiquinol-cytochrome c reductase complex assembly factor 3
 8: ENSG00000072501  0.2871578 7.901466  4.138382 6.440253e-04 0.5074877 -0.34573791    up.reg     8243  SMC1A                   structural maintenance of chromosomes 1A
 9: ENSG00000163346 -0.3361114 6.606192 -4.118766 6.724169e-04 0.5074877 -0.41793193  down.reg    57326 PBXIP1                         PBX homeobox interacting protein 1
10: ENSG00000234444  0.4092223 4.950297  4.180120 5.875733e-04 0.5074877 -0.57590263    up.reg   728927 ZNF736                                    zinc finger protein 736

#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Treatment.arm) + geom_boxplot() + stat_compare_means(aes(label=paste0("p.adj = ",..p.adj..)),method="t.test") + theme_tufte() + theme(legend.position = "top")

## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Treatment.arm)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")

#GSEA - c6
gsea <- GSEA.run(DE.data$v,d,cont.matrix[,1],Hs.c6,0.05)

'select()' returned 1:many mapping between keys and columns

gsea[Direction=="Up",.(rn,NGenes,Direction,FDR)]

                         rn NGenes Direction          FDR
 1:         VEGF_A_UP.V1_DN    167        Up 2.498761e-05
 2:           ERB2_UP.V1_DN    186        Up 4.565218e-05
 3:              TBK1.DF_DN    267        Up 6.914354e-05
 4:           PIGF_UP.V1_UP    156        Up 5.032092e-04
 5:  GCNP_SHH_UP_LATE.V1_UP    157        Up 1.941501e-03
 6:           E2F1_UP.V1_UP    171        Up 4.030628e-03
 7: GCNP_SHH_UP_EARLY.V1_UP    146        Up 4.030628e-03
 8:          RPS14_DN.V1_DN    154        Up 4.522622e-03
 9:           LTE2_UP.V1_DN    177        Up 6.406573e-03
10:      CSR_EARLY_UP.V1_UP    145        Up 7.089887e-03
11:      BCAT_BILD_ET_AL_DN     38        Up 1.263871e-02
12:      PRC2_EZH2_UP.V1_UP    132        Up 1.722478e-02

gsea[Direction=="Down",.(rn,NGenes,Direction,FDR)]

              rn NGenes Direction         FDR
1: ERB2_UP.V1_UP    174      Down 0.007089887
2:      EIF4E_UP     78      Down 0.027986355

Vehicle

rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Treatment.arm=="vehicle"])
treatment.group <- rnaseq.subset$meta[,Group]
d <- model.matrix(~ 0 + treatment.group)
d

   treatment.group1 treatment.group2 treatment.group3
1                 1                0                0
2                 1                0                0
3                 1                0                0
4                 1                0                0
5                 0                0                1
6                 0                0                1
7                 0                0                1
8                 0                0                1
9                 0                0                1
10                0                1                0
11                0                1                0
12                0                1                0
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment.group
[1] "contr.treatment"

cont.matrix <- matrix(c(-1,1,0,-1,0,1,0,-1,1), dim(d)[2], 3)
colnames(cont.matrix) <- c("Long.minus.short","Resistant.minus.short","Resistant.minus.long")
rownames(cont.matrix) <- colnames(d)
cont.matrix

                 Long.minus.short Resistant.minus.short Resistant.minus.long
treatment.group1               -1                    -1                    0
treatment.group2                1                     0                   -1
treatment.group3                0                     1                    1

#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,3,cont.matrix = cont.matrix)

dt <- decideTests(DE.data$efit)
summary(dt)

       Long.minus.short Resistant.minus.short Resistant.minus.long
Down               2243                  2816                    0
NotSig            12186                 11123                16786
Up                 2357                  2847                    0

dt.common <- dt[rowSums(abs(dt))!=0,]
vennDiagram(dt.common)

#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Group) + geom_boxplot() + theme_tufte() + theme(legend.position = "top")

## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Group)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")

Drug

rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Treatment.arm=="GDC0032"])
treatment.group <- rnaseq.subset$meta[,Group]
d <- model.matrix(~ 0 + treatment.group)
d

   treatment.group1 treatment.group2 treatment.group3
1                 1                0                0
2                 1                0                0
3                 1                0                0
4                 1                0                0
5                 0                0                1
6                 0                0                1
7                 0                0                1
8                 0                0                1
9                 0                0                1
10                0                0                1
11                0                0                1
12                0                0                1
13                0                1                0
14                0                1                0
15                0                1                0
16                0                1                0
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment.group
[1] "contr.treatment"

cont.matrix <- matrix(c(-1,1,0,-1,0,1,0,-1,1), dim(d)[2], 3)
colnames(cont.matrix) <- c("Long.minus.short","Resistant.minus.short","Resistant.minus.long")
rownames(cont.matrix) <- colnames(d)
cont.matrix

                 Long.minus.short Resistant.minus.short Resistant.minus.long
treatment.group1               -1                    -1                    0
treatment.group2                1                     0                   -1
treatment.group3                0                     1                    1

#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,4,cont.matrix = cont.matrix)

dt <- decideTests(DE.data$efit)
summary(dt)

       Long.minus.short Resistant.minus.short Resistant.minus.long
Down                 78                  1791                  514
NotSig            16695                 13186                16005
Up                   15                  1811                  269

dt.common <- dt[rowSums(abs(dt))!=0,]
vennDiagram(dt.common)

#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Group) + geom_boxplot() + theme_tufte() + theme(legend.position = "top")

## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Group)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")